function [IA_Network, consistent, iterations] = IA_Network_Path_Consistency_Parallel(IA_Network, i, j, IA_Relation)
%IA_Network_Path_Consistency To Add R(i, j)
%   Detailed explanation goes here

% Todo (RF): both 'N' and 'R' matrices are at least 2x bigger than required 
% as the relationship matrix is symmetric around the (top left to bottom 
% right) diagonal. The diagonal also contains no information as each 
% interval can only be equal to itself and no other relation is possible.

consistent = 1;
iterations = 0;

if IA_Relation_Is_Dummy(IA_Relation)
  return;
end

if IA_Relation_Is_Null(IA_Relation)
  consistent = 0;
  return;
end

N = IA_Network.relations;
N(i, j) = IA_Relation_Intersect(N(i, j), IA_Relation);
N(j, i) = IA_Relation_Invert(N(i, j));

if IA_Relation_Is_Null(N(i, j))
  consistent = 0;
  return;
end

% Add <i,j> to queue ToDo;
ToDo = struct('i', i, 'j', j);

% While ToDo is not empty do
% begin
done = 0;
while (size(ToDo, 2) > 0) && (~done)
  % N^2 = N o N
  newToDo = struct('i', i, 'j', j);
  newToDo = newToDo(2:size(newToDo, 2));
  %parNoN(1:length(ToDo), 1:size(N,1), 1:size(N,2)) = IA_Relation_Dummy();
  for k = 1:length(ToDo)
    curr = ToDo(k);
    i = curr.i;
    j = curr.j;

    % Note: Parallel threads can operate on the same matrix, they must just 
    % be locked out of accessing the same locations at the same time (both 
    % for intervals that are updated as well intervals looked up).
    %
    %[parNoN(k,:,:), consistent, newToDo] = IA_Network_Path_Consistency_Parallel_Iteration(N, i, j, newToDo); %#ok<AGROW>
    [N, consistent, newToDo] = IA_Network_Path_Consistency_Parallel_Iteration(N, i, j, newToDo);
    if ~consistent
      return;
    end
  end
  
% Note: Parallel threads now operate on the same matrix, so the code below
% is no longer required.
%
%   N(:,:) = parNoN(1,:,:);
%   for k = 2:length(ToDo)
%     for i = 1:size(N,1)
%       for j = 1:size(N,2)
%         %N(i, j) = IA_Relation_Intersect(N(i, j), parNoN(k, i, j));
%         N(i, j).set = N(i, j).set & parNoN(k, i, j).set;
%       end
%     end
%   end
%
% Todo: Since I do not have the Matlab parallel toolbox I am only 
% 'simulating' parallel execution. Port to Python and run in parallel.
  
  ToDo = newToDo;
  iterations = iterations + 1;
end

IA_Network.relations = N;
end